% sstate.m  
% establishes parameter values and steady state for the LMV ReSTAT model of oil production and usage 


format short

% PARAMETER VALUES. 

beta = 0.99;                % Discount factor
delta = 0.025;              % Depreciation rate    


psi1 = 0.1056;  %constant controlling n_ss
psi2 = 1;

sigma = 2.0;

hab = 0.0;                  % habit formation if any


% Production structure
omega1 = 0.0585;      % weight of oil in CES composite of VA and oil

omega2 = 0.0585;      % weight of oil in CES composite of goods and oil in Household consumption

omega3 = 0.66;      % weight of labour intput in CES composite of labour and capital to produce VA

eps1 = 0.1;         % elas. of subst. between oil and VA in CES composite for Y production

eps2 = 0.1;         % elas. of subst. between oil and goods in CES composite Household consumption

eps3 = 0.99999;     % elas. of subst. between N and K CES composite for VA production
                                        
oil_ss = 1;  % qty of oil available in ss (scale variable)

% Oil Inventories Sector
kappa = -0.05;   %"convenience yield"
alpha_inv = 0.9999;  % "oil waste"
s_ss  = 0.5;  % s.s. value of inventories as frction of s.s. total quantity available
answinv = param(9);
if answinv == 1
    % economy with storage
    % calibrate cost of inventory holding to rationalize s_ss
    phi_inv = (alpha_inv*beta-1-kappa)/s_ss;
elseif answinv == 0 
    %"no Storage case":
    % Inventories are nil
   s_ss = 0;
   phi_inv = 10000000.0;
end

% capital adjustments costs:
answcost = 1;   %1 for standard, 2 for CEE
if answcost == 1
   phicost = 0.00000001;
   %phicost = 10.1;
elseif answcost ==2
    phi_CEE = 1.0;
end

answfutures = 1;
% set to 1 if you want "model coherent" computations for the futures
% prices. 


trait1 = 'k-';
trait2 = 'k--';
trait3 = 'k:';
epais = 2.0;


% Solving for the steady state value

bigR = 1/beta;
rk = 1/beta-1+delta;

gov = 0.0;

% Find endogenous value of oil's (relative) price, which solves ss
x0 = 2.5;
pofind = fzero('solvess',x0,[],rk,delta,sigma,psi1,psi2,omega1,eps1,omega2,eps2,omega3,eps3,s_ss,phi_inv,alpha_inv,kappa);

if(0)
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') 
disp('                 ')   
disp( [ 'S.S. value of p^oil found using fzero : ' num2str(pofind) ]) 
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') 
disp('                 ') 
end

% With p^oil done, ready to compute rest of stuff
qy = omega1*pofind^(-eps1);
vy = ( 1-omega1^(1/eps1)*qy^((eps1-1)/eps1) ) / (1-omega1)^(1/eps1);
vy = vy^(eps1/(eps1-1));
aide = ((1-omega1)/vy)^(-1/eps1);
vk = (rk*aide*(1-omega3)^(-1/eps3))^eps3;
kn = ( vk^((eps3-1)/eps3) - (1-omega3)^(1/eps3) ) / omega3^(1/eps3);
kn = kn^(eps3/(1-eps3));
vn = vk*kn;
yn = vn/vy;
wage = ((1-omega1)/vy)^(1/eps1) * (vn*omega3)^(1/eps3);
p = (omega2*pofind^(1-eps2)+1-omega2)^(1/(1-eps2));
cn = (yn-delta*kn)/( (1-omega2)*p^(eps2));

aide = p*(cn^(sigma))*psi1/wage;
n = aide^(-1/(sigma+psi2));
    
y = yn*n;
q = qy*y;
c = cn*n;
co = omega2*c*(pofind/p)^(-eps2);
cg = (1-omega2)*c*(1/p)^(-eps2);
po = pofind;
k = kn*n;
va = vn*n;


if(0)
%Display some elements of the s.s.
disp('                 ') 
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')                      
disp(' Steady-State value of individual variables') 
disp('         ')
disp('         ') 
disp( '      k         c          n           y     wages    rk')
disp([k c n y wage rk])

disp('         ') 
disp( '    q^oil    inv^oil  p^oil')
disp([   q s_ss  po])
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') 
disp('                 ')  


disp('                 ')  
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')                     
disp(' Steady-State value of relevant ratios') 
disp([' oil inventories/total supply : ' num2str(s_ss) ])
disp([' Value of oil inventories/GDP : ' num2str(po*s_ss/y)])
disp([' Labour share in output : ' num2str(wage*n/y) ])
disp([' Capital share in output : ' num2str(rk*k/y) ])
disp([' Oil share in output (p^o * q / GDP) : ' num2str(po*q/y) ])
disp([' Oil share in Household Consumption basket (p^o * c^o / p C) : ' num2str(po*co/(p*c)) ])
disp('                 ')                                      
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') 
disp('                 ') 
end


% Setting up parameters related to driving processes
% -------------------------------------------------------

% The (2 by 1) vector of observed exogenous variables is
% x_t = [shock_supply_oil shock_neutral_TFP]'
% The (4 by 1) vector of underlying exogenous states is 
% ksi_t = [shock_supply(pers.) shock_supply(trans.) shock_neutraltechno(persistent) shock_neutraltechno(transitory)]' 
% the system evolves according to
% x_t = H' ksi_t
% ksi_t+1 = F ksi_t + eta_t+1
% which corresponds to KW(2002)'s notation on page 73

nstates = 4;

rhoe_p = param(1);   % serial correlation of persistent energy shocks
sige_p = param(2); % std. deviation of innovations to those shocks
rhoe_t = param(3);   % serial correlation of transitory energy shocks
sige_t = param(4);    % std. deviation of innovations to those shocks

%Calibration for TFP shocks "early 2000s"
rhoz_p = param(5);   % serial correlation of persistent productivity shocks
sigz_p = param(6); % std. deviation of innovations to those shocks
rhoz_t = param(7);   % serial correlation of transitory productivity shocks
sigz_t = param(8);    % std. deviation of innovations to those shocks


 whatinfo = param(10);
 answlearn = param(11);
 


     
F = [rhoe_p  0       0          0       
     0       rhoe_t  0          0       
     0       0       rhoz_p     0       
     0       0       0          rhoz_t  ];

H =  [1 1 0 0   
      0 0 1 1  ]';
 
QQ = [sige_p^2   0            0            0         
         0          sige_t^2     0            0        
         0          0            sigz_p^2     0         
         0          0            0            sigz_t^2  ];
     
     

%starting value for the MSE in the updating ( P(1)|0 )
P = reshape(inv(eye(nstates*nstates)-kron(F,F))*reshape(QQ,nstates*nstates,1),nstates,nstates);


   
   if(1)
   % running this a few times will insure the updating starts with the stationary P 
   % rather than P(1) | 0. It seems that within 15 iterations, we are essentially
   % at the s.s. P
      for i = 1:50
        Ptest = P;
        K = F*P*H/(H'*P*H);             %eq. 13.2.19 Hamilton  
        P = (F-K*H')*P*(F'-H*K')+QQ;        %eq. 13.2.28 Hamilton
      end
   end


   
% For the benefit of working with the program <mdrkw>, redefine matrices
Q = H';
RHO = F;


if whatinfo == 1
    %if information is complete, no need for learning via the
    %Kalman filter: the "gain" is one
    bigA = eye(nstates);
elseif whatinfo ==2
    % when information is incomplete, updating via the Kalman filter
    % "bigA" is what will be sent to mdrkw_learning
    bigA = P*H/(H'*P*H);
    % To build Table 5 we want to use the above parameters to draw the shocks
    % but use a gain from another period. This separates the impacts of changed shocks versus changed learning
    if answlearn == 0
        bigA = gainkeep;
    end
    bigA = bigA*H';
    bigA;
end
       
   
